%title:             insurance model
%author:            neale mahoney
%lastest revision:  6/28/14


%% Setup

clear all; clc;

% cd(<folder directory>);

version = 'v7';

%% Primatives

fprintf('Displaying primatives \n\n');

n_i = 16e3;             %number of individuals
n_c = 1e3;              %number of cost shocks per individual
smooth_type = 'lowess'; %moving, lowess, rlowess
bw = .1;                %bandsmooth for moving average smoother


%% Joint distribution of risk aversion and expected medical costs

%distribution of absolute risk aversion from Table 3 of Handel, Hendel and Whinston (2013)
m = 4.39e-4;
v = (6.63e-5)^2;
mu_alpha = log((m^2)/sqrt(v+m^2));
sigma_alpha = sqrt(log(v/(m^2)+1));		

%gamble interpretation
f = zeros(100,1);
for i = 1:100,
    f(i) = cara_function(0,m,[-100 i],[0 0]);
end

[value index] = min(abs(f));
fprintf('Mean risk adversion makes individual different between 50-50 gamble for (100,%i) and 0 \n', index)

% distribution of health type
m = .9794951;
v = 1.377593^2;
mu_l = log((m^2)/sqrt(v+m^2));
sigma_l = sqrt(log(v/(m^2)+1));

%joint normal
mu = [mu_alpha; mu_l];
rho = 0;			%assume risk aversion and health risk are uncorrelated in the population
cov = rho*sigma_alpha*sigma_l;
sigma = [sigma_alpha^2 cov; cov sigma_l^2];
X = mvnrnd(mu,sigma,n_i);

%distribution of risk preferences and health type
alpha = exp(X(:,1));
lambda = exp(X(:,2));

%distribution of realized costs
m = 3138.765; 
v = 10125.74^2;
mu_c = log((m^2)/sqrt(v+m^2));
sigma_c = sqrt(log(v/(m^2)+1));
rho = .498481; 

%distribution of realized cost conditional on health type
c = zeros(n_i,n_c);
part_1 = mu_c;	
part_2 = rho*sigma_c*(X(:,2)-mu_l)/sigma_l;	

for i = 1:n_i   
    part_3 = sqrt((1-rho^2))*sigma_c*randn(1,n_c);
    c(i,:) = exp(part_1 + part_2(i) + part_3);
end

fprintf('Mean absolute risk aversion parameter is %f \n',mean(alpha));
fprintf('Std dev of absolute risk aversion parameter is %f \n',std(alpha));

fprintf('Mean population cost is %i \n',round(mean(reshape(c,n_i*n_c,1))))
fprintf('Std dev of cost is %i \n',round(std(reshape(c,n_i*n_c,1))))

%av 60 plan: low quality plan that covers 60% of medical cost on average
c_oop_60 = min(min(c,500+.4*(c-500)),8e3);  %500 deductible, 8000 out of pocket maximum
c_ins_60 = c - c_oop_60;		
c_ins_mean_60 = mean(c_ins_60,2);	

fprintf('AV 60 mean plan cost is %i \n',round(mean(c_ins_mean_60)))
fprintf('AV 60 mean oop cost is %i \n',round(mean(reshape(c_oop_60,n_i*n_c,1))))
fprintf('AV 60 AV is %f \n',mean(c_ins_mean_60)/mean(reshape(c,n_i*n_c,1)))

% av 90 plan: high quality plan that covers 90% of medical cost on average
c_oop_90 = min(.1*c,8e3);   %no deductible, 8000 out of pocket maximum
c_ins_90 = c - c_oop_90;
c_ins_mean_90 = mean(c_ins_90,2);

fprintf('AV 90 mean plan cost is %i \n',round(mean(c_ins_mean_90)))
fprintf('AV 90 mean oop cost is %i \n',round(mean(reshape(c_oop_90,n_i*n_c,1))))
fprintf('AV 90 AV is %f \n',mean(c_ins_mean_90)/mean(reshape(c,n_i*n_c,1)))

%allow bankruptcy to cap out of pocket risk
c = min(c,20e3);

%% wtp
fprintf('\n\nEstimating WTP \n\n');

v_grid = 0:10:5e4;
v = zeros(n_i,1);
for i = 1:n_i,
    %v: willingness to pay additional premium such that 
    %the customer is indifferent between high and low quality plan
    if mod(i,25) == 0, fprintf('Estimating WTP for observation %i \n',i); end
    v(i) = wtp_function(alpha(i),v_grid,c_oop_90(i,:),c_oop_60(i,:));
	
end


%% Equilibrium - Baseline

%sort matrix X in descending order by wtp and smooth
wtp_shift = 900;
X = [v + wtp_shift c_ins_mean_90-c_ins_mean_60];
X = sortrows(X,-1);

v_sort = smooth(X(:,1),bw,smooth_type);	
mc = smooth(X(:,2),bw,smooth_type);
q = transpose((1:n_i)./n_i);

plot_switch = 1;
plot_number = 1;
plot_dir = '';
plot_name = 'all_theta';	

%market power: theta = 0.5
%this corresponds to Cournot competition with two high-quality plans
theta = .5;
out_theta = equilibrium_function_v2(...
    v_sort,mc,q,theta,...
    bw,smooth_type,plot_switch,plot_number,plot_dir,plot_name,version); 

%% Equilibrium - Risk adjustment
%sort matrix X in descending order by wtp and smooth
wtp_shift = 900;
X = [v+wtp_shift c_ins_mean_90-c_ins_mean_60];
X = sortrows(X,-1);

v_sort = smooth(X(:,1),bw,smooth_type);
mc = smooth(X(:,2),bw,smooth_type);
q = transpose((1:n_i)./n_i);

sigma = 0;	%perfect risk adjustment
ac_bar = mean(mc(:,1));	
mc_ra = sigma*mc + (1-sigma)*ac_bar;	

plot_switch = 1;
plot_number = 2;
plot_dir = '';
plot_name = 'ra_theta';

theta = .5;
out_theta = equilibrium_function_v2(...
    v_sort,mc_ra,q,theta,...
    bw,smooth_type,plot_switch,plot_number,plot_dir,plot_name,version); 

%% Equilibrium - Negative Risk adjustment

%sort matrix X in descending order by wtp and smooth
wtp_shift = 900;
X = [v+wtp_shift c_ins_mean_90-c_ins_mean_60];
X = sortrows(X,-1);

v_sort = smooth(X(:,1),bw,smooth_type);
mc = smooth(X(:,2),bw,smooth_type);
q = transpose((1:n_i)./n_i);

sigma = 2;  %negative risk adjustment
ac_bar = mean(mc(:,1));
mc_ra = sigma*mc + (1-sigma)*ac_bar;	

plot_switch = 1;
plot_number = 3;
plot_dir = '';
plot_name = 'neg_ra_theta';

theta = .5;
out_theta = equilibrium_function_v2(...
    v_sort,mc_ra,q,theta,...
    bw,smooth_type,plot_switch,plot_number,plot_dir,plot_name,version); 

%% Equilibrium -  Segmented Markets
 
fprintf('\n\nEstimating equilibrium in segmented markets \n\n');

r = lambda;	%health type

n_segments = 4;
bw = .2;

%split the market by health type
if n_segments == 2,
    cut_points = median(r);
else 
    cut_points = quantile(r,n_segments-1);
end

cut_points = [min(r) cut_points max(r)];
segment = struct([]);

for i = 1:n_segments,
    
    X = zeros(1,2);
    index = 0;
    for j = 1:length(r),
        if r(j) > cut_points(i)  && r(j) <= cut_points(i+1)
            
            index = index + 1;
            wtp_shift = 900;
            X(index,:) = [v(j) + wtp_shift c_ins_mean_90(j)-c_ins_mean_60(j)];
        end
    end
    
    %sort matrix X in descending order by wtp and smooth
    segment(i).X        = sortrows(X,-1);
    segment(i).v_sort   = smooth(segment(i).X(:,1),bw,smooth_type);
    segment(i).mc       = smooth(X(:,2),bw,smooth_type);
    segment(i).q        = transpose((1:index)./index);
    
    fprintf('Segment %i runs from %f to %f. It has %i obs and mean cost of %f \n',...
        [i cut_points(i) cut_points(i+1) index mean(segment(i).mc)]); 
end

%plot graphs for each segment, with market power theta = 0.5 no adjustment of selection
plot_switch = 1;
plot_dir = '';
out_segment = zeros(8,n_segments);
n_i_segment = zeros(n_segments);

for i = 1:n_segments,
    
    plot_number = 4 + i;
    plot_name = ['segment_' num2str(i)];
           
    theta = .5;
    out = equilibrium_function_v2(...
        segment(i).v_sort,segment(i).mc,segment(i).q,theta,...
        bw,smooth_type,plot_switch,plot_number,plot_dir,plot_name,version);
    
    n = length(segment(i).q);
    
    out_segment(:,i) = [out(:,2); n]; 
end

